Influence of vortices and phase fluctuations on thermoelectric transport properties of 

superconductors in a magnetic fleld 
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We study heat transport and thermoelectric effects in two-dimensional superconductors in a mag- 
netic field. These are modeled as granular Josephson-junction arrays, forming either regular or 
random lattices. We employ two different models for the dynamics, relaxational model-A dynamics 
or resistively and capacitively shunted Josephson junction (RCSJ) dynamics. We derive expressions 
for the heat current in these models, which are then used in numerical simulations to calculate 
the heat conductivity and the Nernst coefficient for different temperatures and magnetic fields. At 
low temperatures and zero magnetic field the heat conductivity in the RCSJ model is calculated 
analytically from a spin wave approximation, and is seen to have an anomalous logarithmic depen- 
dence on the system size, and also to diverge in the completely overdamped limit C — >■ 0. From our 
simulations we find at low magnetic fields that the Nernst signal displays a characteristic "tilted 
hill" profile similar to experiments and a nonmonotonic temperature dependence of the heat con- 
ductivity. We also investigate the effects of granularity and randomness, which become important 
for higher magnetic fields. In this regime geometric frustration strongly influences the results in 
both regular and random systems and leads to highly nontrivial magnetic field dependencies of the 
studied transport coefficients. 

PACS numbers: 74.81.-g, 74.25.F-, 74.25.fg, 74.25.Uv, 74.78.-w 



I. INTRODUCTION 

Thermoelectric effects in superconductors are of con- 
siderable interest, since they provide an important probe 
of fluctuations and correlations in these materials. Such 
effects have gained an increasing amount of attention 
since the recent measurements of the Nernst effect in the 
pseudogap phase of underdoped high- Tc cuprates, where 
a remarkably large effect was observed. ^ The Nernst ef- 
fect has since then been measured in many other ma- 
terials, e.g., in superconducting thin films^ and in the 
iron-pnictideSiS: The Nernst effect is usually very small 
for conventional metals, making it a particularly sensitive 
probe of superconducting fluctuations. Theoretical and 
numerical studies have described the large Nernst effect 
either in terms of superconducting fluctuations of Gaus- 
sian nature;^"" or as fluctuations of the phase of the order 
parameter, i.e., as a result of vortices J^i^ Other expla- 
nations of nonsuperconducting origin have also been put 
forward, e.g., proximity to a quantum critical point /'^i^^ 
quasiparticlesr^ii^ stripe order^ etc. Here we will fo- 
cus exclusively on the effects of phase fluctuations and 
vortices on heat and charge transport. Phase fluctua- 
tions were early on proposed to play a key role in the 
pseudogap phase of underdoped cuprates In quasi- 
two-dimensional superconducting films and Josephson- 
junction arrays they are known to be the dominant fluc- 
tuations.^- 

In this paper we study the heat transport and the 
thermoelectric response in two-dimensional granular su- 
perconductors or Josephson- junction arrays, using two 
widely employed models for the dynamics, (i) relax- 
ational model-A dynamics, described by a Langevin 
equation, or (ii) resistively and capacitively shunted 



Josephson junction (RCSJ) dynamics. These are well 
suited for numerical simulation studies, and have been 
used extensively to calculate electric and magnetic static 
and dynamic properties, and to study vortex dynam- 
ics under influence of electric currents Ji"— The cal- 
culation of thermoelectric properties is, however, less 
straightforward. Previous simulation studies have used 
time-dependent Ginzburg-Landau theory,— phase-only 
XY models^ with Langevin dynamics, or vortex dynam- 
icsi^ They have been limited to rather narrow parameter 
regimes with a focus on explaining the large Nernst ef- 
fect seen in the pseudogap phase of underdoped cuprates. 
Here we present a comprehensive study of heat conduc- 
tivity, thermoelectric effects, and resistivity for a broad 
range of parameters. We also investigate the effects of 
a granular structure. The models are defined on a dis- 
crete lattice and can be experimentally realized in granu- 
lar superconductors and artificially fabricated Josephson- 
junction arrays. At low magnetic fields discreteness ef- 
fects become less important so that our results in this 
regime are relevant also for homogenous bulk supercon- 
ductors. On the other hand, in granular superconductors 
transport properties are strongly affected by discreteness 
and geometric frustration effects. This is particularly 
true at high magnetic fields when the inter- vortex dis- 
tance becomes comparable to the granularity. This leads 
to a rich structure in, e.g., the Nernst signal as the mag- 
netic field is varied, with anomalous sign changes occur- 
ring in the vicinity of special commensurate fillingsi^ 



To start with, let us first recall that the heat current 
density J'^ and electric current density J are related to 
the thermodynamic forces, the electric field E and the 
temperature gradient — VT, by the standard phenomeno- 
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logical linear relations 
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where the thermoelectric and electrothermal tensors obey 
the Onsager relation a = Ta. The Nernst coefficient v is 
defined as the ofi'-diagonal response of the electric field Ey 
to an applied temperature gradient V^jT in a transverse 
magnetic field B^, 
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where is the so called Nernst signal. Both the Nernst 
effect and the heat conductivity are measured under the 
condition of zero electric current, such that 



Cat 



■*-xx^ xy 



aa a. 



(3) 
(4) 



In a system with no Hall effect {axy = 0), which is the 
case for the models studied below, Eq. ^ reduces to 
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In a phase-fluctuating superconductor, mobile vortices, 
either thermally excited or induced by an applied mag- 
netic field, may significantly affect transport properties. 
An applied electric current J will exert a perpendicular 
force on the vortices and their motion will generate a 
transverse electric field E = — v x B parallel to J, thus 
causing a finite flux-flow resistivity. Vortex motion can 
also be caused by a temperature gradient, thereby gen- 
erating an electric field perpendicular to both the mag- 
netic field and the temperature gradient. The Nernst 
coefficient defined in Eq. © can be seen simply as the 
diagonal response v = —vVT of the vortex velocity to 
the temperature gradient. For this reason it is plausible 
that a large Nernst signal is expected in the vortex liquid 
phase. 

Via an Onsager relation a heat current can then be 
generated from an applied electric current. The vor- 
tices therefore also contribute to the heat conductivity, 
although other heat carriers, e.g., phonons or quasiparti- 
cles, probably dominate. From the vortex point of view it 
is natural to consider the applied current as the driving 
force and the electric field, which is proportional (but 
perpendicular) to the vortex current, as the response. 
This yields an alternative, but equivalent formulation of 
the linear relations in Eq. ([T]) 
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where /? = —P/T — —a^^a. This is the approach em- 
ployed in our simulations. Instead of calculating cn 
through other transport coefficients, i.e., using Eq. ([3]), 
we obtain the Nernst signal directly as cn — Pyx — 
Pxy/T for J = 0. The longitudinal heat conductivity is 



just the diagonal components of the tensor k in Eq. ([5|), 
and p = (T~^ is the resistivity. 

The picture described above applies when the vortices 
are free to move in response to the driving forces. Pinning 
of vortices to material defects, grain boundaries, etc., can 
lead to dramatic changes of the transport properties. 

The paper is organized as follows. In Sec. |ll]we intro- 
duce the models and their dynamics, Langevin or RCSJ, 
on general two-dimensional (2D) lattices. In Sec. IIIII we 
derive an expression for the heat current for the models 
studied. This has over the years proven to be a subtle 
task, especially in the presence of a magnetic field. We 
present two separate routes to finding the explicit form of 
the heat current for Langevin and RCSJ dynamics. In ad- 
dition we show, using a functional integral approach, that 
the Nernst signal indeed can be calculated from a Kubo 
formula involving the heat current. Section ITVl discusses 
the thermal conductivity at zero magnetic field in the low 
temperature limit, where a spin wave approximation is 
applicable. Our analytic calculations reveal a logarithmic 
system size dependence of k in this regime for the RCSJ 
case. In Sec. |V]we give some technical details of our nu- 
merical simulations, and in the last part. Sec. I VII the 
results of our numerical simulations on square and irreg- 
ular lattices are presented. We consider the case of zero 
and weak magnetic fields as well as the high field limit, 
where the transport coefficients are severely affected by 
geometric frustration. In the weak magnetic field limit 
the results are discussed in relation to previous theoreti- 
cal works and experiments. 



II. MODELS 

We model a 2D granular superconductor (of size LxL) 
as an array of superconducting grains connected by 
Josephson junctions. These grains may or may not be 
ordered and connected in a regular fashion. The super- 
current flowing between two superconducting grains i and 
j is given by the Josephson current-phase relation 
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(6) 

A • dr, (7) 



where If^ is the critical current of the junction, $o = h/2e 
is the superconducting flux quantum, and 9i is the super- 
conducting phase of grain i. The grains are assumed to 
be small enough (< the coherence length) that the phase 
is constant over each grain. Further, A is the vector 
potential, which we here separate into two parts 



A(r,i) = Aext(r) + ^A(t), 



(8) 



where Aoxt(r) is constant in time and corresponds to a 
uniform magnetic fleld B = V x A perpendicular to the 
array, and A(i) = (Aj.(i), Aj,(t)) is spatially uniform, 
but time dependent and is necessary to include in order 
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to describe temporal fluctuations in the average electric noise in the resistors Rq, have the correlations 



field E = — ^A, when periodic boundary conditions are 
usedi ^^i^" Local fluctuations in the magnetic field B and 
hence A are ignored. 



iv^it)) = 0, {v^{tH{t')) - ^^S,,S{t - t'), (14) 

(C(<)> = 0, {Ut)Ut')) = ^^^^P^s^Jit - 1'). (15) 



A. Langevin dynamics 

The Langevin dynamics represents a phase-only de- 
scription of the time-dependent Ginzburg-Landau dy- 
namics (TDGL) in uniform magnetic field. The phe- 
nomenological equations of motion for {9} and A are 
of local relaxation type 



H = 



where the time constant 7 is dimensionless and 7a = 
7L^, H is the XY model Hamiltonian, and rji and C are 
Gaussian stochastic white noise processes. 

These equations can be cast in the form of circuit equa- 
tions for an electric circuit built up using the elements 
displayed in Fig. [U where each site i is connected via a 
resistor Rq — h/Ae^^ to ground: 
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The sum in the first equation is taken over the set Mi 
of superconducting grains connected to i and is equiva- 
lent to a lattice divergence. In the second equation {ij) 
denotes a sum over all links in the system. The vector 
Yji = Yj Yi goes from site i to site j and J"'^* is a fixed ex- 
ternal current density. The Gaussian white noise terms rji 
and ^, which now can be interpreted as Johnson-Nyquist 



B. RCSJ dynamics 



In the RCSJ model each Josephson junction with crit- 
ical current If^ is shunted by both a resistor Rij and a 
capacitor Cij, see Fig. [51 The RSJ model is obtained as 
a special case when setting Cij = 0. We write the to- 
tal current from site i to site j as a sum of the super-, 
resistive, capacitative, and noise currents 
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where the voltage between grain i and j is 

■ h 

V^j = 14 - Vj - Ay = — jij, 



(16) 



(17) 



and the Johnson-Nyquist noise in the resistors has zero 
mean (/j" ) = and covariance 



(J2W4"z(i')> 



2kBT 

Rii 



{5,k5ji~5u5,kW~t')- (18) 



The equations of motion for the phases {9i} and the 
twists A are obtained from local current conservation at 
each grain 



(19) 



Here we introduced a "transport current" which is 
a static, divergence free current distribution, satisfying 
any external boundary conditions, but is otherwise ar- 
bitrary. One may, for instance, connect some nodes to 
fixed external current sources or sinks. In the present 
work we will usually use periodic boundary conditions 
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FIG. 1. Equivalent circuit model for Langevin dynamics. The 
cross denotes the Josephson junction with critical current I^j. 



FIG. 2. Equivalent circuit model for RCSJ dynamics. 
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instead, together with the condition of a fixed average 
current density J°'^' in the system, 

^/^'r,,^^/^r,,=L2jcxt^ (20) 

For model purposes we may define a local current density 
on the links of the lattice as 

JW=E r nT5{r-^')dv'. (2f) 

which directly leads to Eq. pO| when averaged over the 
system. 




FIG. 3. (Color online) A part of a random lattice with the 
length dfj and the vector Vji defined. The green lines make 
up the direct Delaunay lattice and the orange structure is the 
dual Voronoi lattice. 



C. Lattice structure 

We are interested in modeling both regular and ran- 
dom granular superconductors. At low magnetic fields 
the vortex separation is large compared to the granules 
and the lattice structure does not matter much. In this 
regime the models approximate bulk superconductors. 
In the opposite limit of high magnetic fields the lattice 
structure is important as frustration effects strongly in- 
fluence the transport properties. Note that the formula- 
tion of the models defined above is independent of lat- 
tice structure. We will limit ourselves to two dimensions 
in the present work. In our simulations presented be- 
low in Sec. I VII we consider not only square lattices, but 
also random lattices appropriate as models of random 
granular superconductors. These irregular lattices are 
constructed by generating a set of randomly distributed 
points Vi = (xiyUi) with unit density on a square and 
connecting nearest neighbors by Delaunay triangulation. 
To control the randomness we use the parameter dmim 
which is the shortest allowed distance between any two 
points in the system. Large values of c?min will thus cre- 
ate a relatively ordered lattice structure, whereas lattices 
with small dmin values are more disordered. For example, 
a given value of dmin = {0.0,0.4,0.6,0.8} corresponds to 
a distribution of grain size diameters with standard devi- 
ation {0.30, 0.23, 0.15, 0.08}, respectively. Examples can 
be seen in Fig. [Tsl For the random lattices we use two 
different models, one where the critical current of every 
junction is set to a constant Ifj = I'^ and one where the 
critical current is proportional to the contact length dfj 
between the grains, I^^ ^ djj (in the RCSJ case we also 
take ~ djj and Cij ^ d-^). This length is simply the 
distance between the points in the dual Voronoi lattice, 
see Fig. H 

D. Transport coefflcients 

In the models defined above the only nonzero transport 
coefficients are the Nernst signal ejv, the diagonal ther- 
mal conductivity k and the electrical resistivity p. These 



may be obtained from equilibrium correlation functions 
using Kubo formulas i ^^i^^ While p relates E to a me- 
chanical perturbation, cm and k give the response to a 
nonmechanical property, namely a temperature gradient, 
and the applicability of a standard Kubo formula is not 
evident. Nevertheless, we show in the next section that 
the transport coefficients can be expressed in standard 
form as 

Q r°° 

ew = ^y {Ey{t)J^{Q))dt, (22) 

^^k^l {J?it)J?mdt, (23) 

P^k^l {Ey{t)EyiO))dt, (24) 

where = LxLy is the area of the system. In these 
equations Ey = — ($o/27r)Ay is the average electric field 
in the y-direction and 

•^-^ = ?) E (^M^^ + ^o) - 4-4) - n^) (25) 

is the average heat current in the x-direction, with Xji — 
Xj — Xi and xlj — {xi + Xj)/2. This form of the heat 
current, which is one of our main results, will be explicitly 
derived in the next section for the specific models under 
consideration. 



III. HEAT CURRENT 

In order to calculate thermoelectric effects and the heat 
conductivity an expression for the heat current is needed. 
Several microscopic derivations of the heat current in su- 
perconductors has been given in the literature. 2^"38 xhe 
presence of a magnetic field yields a complication in that 
the total energy and charge currents consist of magneti- 
zation currents in addition to the transport currentsi ^^i'^° 

Rather than relying on these microscopic expressions 
we derive here, within the framework of the models de- 
fined in Sec. |lTl an expression for the heat current that 
can be used consistently in our calculations. Below we 
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will first derive the energy current by considering the 
continuity equation for the energy density. This gives 
the heat current after subtracting the nondissipative en- 
ergy transport, except for a magnetization contribution, 
which is calculated separately in Sec. IIII A31 

Since we are interested in the response of the system to 
an applied temperature gradient, which is a nonmechan- 
ical thermodynamic force, the standard derivation of a 
Kubo formula does not hold. One approach is to follow 
Luttinger;^ and introduce a "gravitational" field, which 
couples to the energy density, and then proceed with the 
calculation of the response to perturbations in that field. 
For the present models, however, there is an alternative 
route. With the dynamics governed by stochastic equa- 
tions, Eqs. (|12m3p or ()19II20|) . in which the temperature 
enters via the strength of the stochastic noise, it is pos- 
sible to introduce local temperature variations, such as 
a temperature gradient, and calculate the resulting re- 
sponse. This calculation, done in Sec. IIII Bl gives us both 
the Kubo formula Eq. ([^ and the heat current Eq. (pSj) . 
Note that in this setting, each point in the model, or more 
precisely, each resistor in the circuit, is in contact with 
a local heat bath. For finite gradients one cannot expect 
the heat current to be automatically conserved, since the 
resistors act as sinks and sources. For an infinitesimal 
thermal gradient the heat current will, however, be con- 
served on average. Alternatively, one could adjust the lo- 
cal temperatures to make sure that heat is conserved on 
average also for finite temperature gradients. Such self- 
consistent temperatures have been employed in studies of 
heat conductivity in harmonic latticesi^ For the problem 
at hand, however, the temperature profile is determined 
by the total heat transport, including phonons, etc, so 
an externally imposed temperature gradient is probably 
more realistic. For the linear response the form of the 
profile should not matter as long as it is smooth. We will 
use a linear temperature profile below. 

A. Heat current from continuity equations 

In this section we will derive the heat current expres- 
sions for granular superconductors described by Langevin 
and RCSJ dynamics, as used in our simulations. First, 
however, it is useful to discuss briefly the continuum for- 
mulation. Starting from the thermodynamic relation 

Tds = de- fidn - H • dB - E • dD (26) 

for a superconductor in a magnetic field H and electric 
field E, where s and e are the entropy and energy densi- 
ties, fj, the chemical potential, and n the density of charge 
carriers (with charge q), one obtains for the heat current 
density 

JQ = Jg, - - E X H, (27) 

q 

where J is the electric transport current density. The 
total energy current density is the sum of two parts, a 



nonmagnetic part and Poynting's vector, 

Jf„, = + -E X B, (28) 
Mo 

and the transport heat current density therefore also has 
two contribution s "^'^ i ^^ '*^^ 

= - + E X M, (29) 

where M = B//io — H is the magnetization. The latter 
part can be rewritten as 

Ex M = {-V(l)-A) X M 

= -W x(j>M + 0Jm - a X M, (30) 

where J^v/ = Jtot — J = V x M is the magnetization cur- 
rent, and 4> the electric potential. The first term on the 
second line is purely rotational and will not contribute to 
the heat transport when integrated over the system. The 
second term (/'(Jtot — J) combined with the nonelectro- 
magnetic part is the standard expression — {S^/q)3, 
where ^ = /i -I- is the electrochemical potential. The 
last piece — A x M is with our gauge choice [Eq. dH)] spa- 
tially constant. Therefore it cannot be uniquely deter- 
mined via the continuity equations below, but will have 
to be added separately later. 

Note the dual role played by the subtraction E x H 
in Eq. (|27|) . It contains the subtraction <j)J (but in a 
gauge invariant way). At the same time it subtracts the 
magnetic field energy transported by the vortex current, 
since i?z$o can be viewed as a magnetic contribution 
to the vortex chemical potential, while the electric field 
E is proportional (but transverse) to the vortex current 
J„ = E X z/$o- In principle one should also subtract 
the nonelectromagnetic energy transported by vortices 
IJ,v3v In the present models, however, Hy = 0. In fact, 
the chemical potential fi for the charge carriers is also 
zero. We will now derive the analogous expressions for 
the discrete models. 



1. Langevin dynamics 

It is instructive to study a slightly more general model 
which includes the charging energy of the grains, de- 
scribed by a circuit model with capacitors Co added in 
parallel with the resistors Rq to ground. This modifica- 
tion, besides being more general, makes the derivation 
more physically transparent, while leaving the final re- 
sults unchanged. The total energy for this model can be 
written 

^ = - E '^o^ ^'j- + E (31) 

with = 9i — 6j — ^Aij being the gauge invariant 
phase difference between site i and j and Vi = h9i/2e the 
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voltage at site i to ground. With the site i we associate 
a local energy 



The time derivative of this is 



sin7y + CqV^V^. 



(32) 



(33) 



In the first term we may identify Jij'jij sin7ij = Vijlfj. 
The last term contains the current through the capacitor 



tCq 



CoVi, which is eliminated using Kirchhoff's law 



(34) 



where J^^o+" is the current through the resistor to 
ground, including the noise current. This results in 



Ro+n,\ 



ieA/'i 



(35) 



This is on the form of a continuity equation for the local 
energy, since the sum over neighbors j € Afi is the lattice 
analogue of a divergence, and the source term Vil^"'^"' 
represents the dissipated work done by the system on the 
environment. The term involving the vector potential Aij 
will be cancelled when the magnetization contribution is 
deah with in Sec. lIIIA3l The energy current is identified 
as 



(36) 



and by subtracting the transport current we obtain the 
heat current 



(37) 



for Langevin dynamics (excluding the magnetization con- 
tribution) . 



2. RCSJ dynamics 

As in the Langevin case it is convenient to add a ca- 
pacitance Co to ground to the usual RCSJ model. The 
energy for such a model is 

E=-J2 J'^'j cos + + E ^^0^^'' (38) 



where Vij = h"fij/2e = Vi — Vj— Aij is the voltage across 
the junction between site i and j. This implies a local 
energy of the form 

= ^ E ^o^^'J- + + V^.', (39) 

with a time derivative 

= i ^ (jy7ij sin7y + CijVijVij^ + CoViVi. (40) 

The last term is, as in the Langevin case, eliminated using 
current conservation 



E(^+^ 



Il-) + ir = 0. (41) 



where the supercurrent Jf^ = Jf^ sin 7^ , the current 
through the shunting resistor (including the noise) is 
/j^^", the parallel capacitance current if^ = CijVij, and 
the current through the capacitance to ground I^" = 
CoV,. We get 



(42) 



and by introducing the total current — Ifj + I^^"^ 
I^j flowing from i to j, and rearranging 



ei+ 



r 2 



jR+n 

(43) 

These terms have interpretations similar to the Langevin 
case. The second term on the left hand side is the lattice 
divergence of the energy current, the first term on the 
right hand side will be cancelled by the magnetization 
contribution, and the last one represents the work done 
on the environment. The energy current is thus 



and the heat current 



(44) 



(45) 



for RCSJ dynamics (again excluding the magnetization 
contribution). The result is very similar to the Langevin 
case, except that the total current appears instead of just 
the supercurrent. 



3. Magnetization contribution to the heat current 

The models defined in Sec. are formulated in the 
limit where fluctuations in the vector potential are com- 
pletely suppressed, except for the uniform part ^ A. 
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Even so, the latter will, perhaps surprisingly, contribute 
to the heat current. We will derive this contribution in a 
more general setting where local fluctuations in the vec- 
tor potential A are allowed. To write down the magnetic 
energy, we first split the total current flowing between 
sites i and j into transverse and longitudinal parts 

IT ^ IT, + f^i - f^J + A. - A, (46) 

(this is the lattice analogue of writing a vector field as 
the sum of a gradient and a curl). The variables /x/ are 
defined on the plaquettes of the lattice, i.e., on the dual 
lattice, and are often referred to as loop currents, see 
Fig. |4l The remaining part Xi — Xj is the loop current 
of the loop i ^ j ground i, which can be nonzero 
in the presence of resistors Rq and/or capacitors Co to 
ground. Without loss of generality we set Ifj = 0, so that 
the heat and energy currents coincide. The loop currents 
can be used to express the magnetic fluxes through the 
corresponding loops, via the self- and mutual inductances 
of the equivalent circuit diagrams (Fig. [T]-[2]), 

Je^fI J 
^ijo = = C,j (Xi-Xj), (48) 

where we have chosen our gauge such that Ako = 0, where 
denotes the ground, and $™* is the applied flux. The 
sum over J G Afj is a sum over adjacent plaquettes sep- 
arated by the link (ij), with i and j defined by Fig. ID 
With this the magnetization energy is 

To simplify the argument we only included the self in- 
ductances Cij of the loops connecting the ground in the 
last term. 

We may associate a local energy e/ with the plaquette 
/, whose time derivative is 

ei = + \ ^^io(Ai - Xj) 

= ^ E (w - Mj + Aj - Xj) + Aijim + Mj) 
= E E \A^,{^I + ^^J)■ (50) 

The first term in the last line cancels exactly the corre- 
sponding contribution in Eqs. ([35]) and (|43|) . The sum- 
mation in the last term represents the energy flowing into 
the plaquette / from the adjacent plaquettes, and we can 
therefore identify the magnetization contribution to the 
energy current and hence the heat current 

^/^j = -^^(m/ + Mj)- (51) 



It is not hard to see that this result holds also in presence 
of an electric transport current, if the /l(/:s are defined as 
in Eq. (|46l) . Equation (fSTj) is the discrete analogue of the 
last term of Eq. ([50]) . 



^. Full heat current 

Averaging Eq. ([37| or (|45|) over the system and adding 
the magnetization contribution Eq. (j5ip gives the full 
average heat current density in the x-direction 

4^=^E^^4(^^ + ^^^^^^'"^^^ 

E [(4 - ^-r) " «j - ^^j) ^'A ^y"' (52) 

where the x-component of the difference 

vector from site i to j, xlj — {xi-\- Xj)/2, and xj the pla- 
quette centers (the vertices of the Voronoi graph). In the 
limit where A is spatially uniform the last term simpli- 
fies to —AyMz, where Mz = fijili/il is the average 
magnetization density [Vlj is the area of the (Voronoi) 
cell /], in agreement with Eq. pop . A more practically 
useful form is obtained below in Scc. lIII Er4l bv expressing 
the loop currents fii in terms of the total currents /*°*, 
which results in the expression Eq. ([25]) . 



B. Heat current from Kubo formula 

It is useful to see how the heat current enters the lin- 
ear response arising from an applied temperature gradi- 
ent via a Kubo formula. We will do this starting from 
the stochastic equations of motion of either Langevin 
[Eq. (dll)] or RCSJ [Eq. (HH)] dynamics. The temper- 
ature enters these equations only via the strength of the 
Gaussian noise correlation function, which now gets a 
spatial dependence T{x) — Tq — T'xM- T' is considered 




FIG. 4. (Color online) On each plaquette I and J, defined 
with respect to the direct lattice points i and j, there are 
loop currents ^/ and The difference between these loop 
currents ^ij = jji — fij is the transverse part of the total 
current on the link from i to j. 
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a weak perturbation, and we are interested in calculat- 
ing the response of some dynamical variable A to such 
a perturbation. We assume open boundary conditions 
in the x-direction, i.e., in the direction of the applied 
temperature gradient, while in the transverse y-direction 
they can be arbitrary. With no net electric current flow- 
ing through the sample the heat current is equal to the 
energy current. 

For concreteness we consider the response of the elec- 
tric field Ey perpendicular to the temperature gradient 
and the magnetic field, which gives the Nernst signal 
gn = S (Ey) /6T'. Our goal is to express this linear re- 
sponse using a Kubo formula 



T2 



{Eyit)jQ{0))dt, 



(53) 



where then can be identified with the heat current 
(the overbar denotes a spatial average), and D, = L^Ly 
is the system size. 

To this end is convenient to reformulate the problem 
using functional integrals'^'* and write ensemble averages 
as 

(A(t)) = i I Aieit))e-^^'lm[d9], (54) 

where is the statistical weight of a given realization of 
the stochastic process 9{t). The Jacobian J[9] in Eq. ((54)) 
comes from the transformation from the Gaussian noise 
C to 6',<2i but since it plays no role in the following it 
will be dropped henceforth. The linear response of A to 
a temperature gradient can now be expressed (assuming 
time translation invariance) as 



where 



R{t - t') 



Q{t') 



HA{t)) 
ST'{t') 



{A{t)Q{t')) 



SS[0] 



ST'{t') 



(55) 



(56) 



T'=0 



We are interested in the response to a static perturbation 
in a stationary state given by 



R = 



{A{t)Q{0))dt. 



(57) 



Despite the similarity with Eq. (j53l) . we cannot imme- 
diately identify 2T'^Q{t)/n with J'^{t) at this stage, be- 
cause of their different symmetry properties, an issue we 
now digress into. 



1. Symmetry relations of transport coefficients and 
correlation functions 

Before continuing the derivation of the heat current 
from the Kubo formula we discuss some general prop- 
erties of linear response in classical statistical systems. 



It is well known that time reversal symmetry implies 
symmetry relations among the transport coefficients and 
among correlation functions. In the presence of a mag- 
netic field H the time reversal operation should re- 
verse also that. Consider now the correlation function 
CE,Q{t,H) = {Ey{t)Q{0)) entering Eq. which in 

general also depends on the magnetic field H. The elec- 
tric field Ey is even under time reversal, while Q is neither 
even nor odd. Instead we may split Q = Qe + Qo into 
even and odd parts. By parity symmetry the correlation 
functions CsyQ^ „ are odd in H, so that 

Ce,q, {t, H) - Ce,q^ i~t, -i/) - ^Ce,q^ i~t, H), (58) 
CEMi'H) - -CE,QA-t,-H) = CE,QA-t,H). (59) 
This leads to the following symmetry of the response 

R{t, H) = Ce,q^ {t, H) + Ce,q„ {t, H) 

= -CE,QA-t,H) + CE,QA-t,H). (60) 

For t < causality implies that R{t, H) — 0, hence 
CEyQo{—t, H) — CEyQ^{—t, H), so that the response can 
be written solely in terms of the odd part of Q as 



R{t,H)=2e{t)CEyQAt^H), 



(61) 



where Q{t) is the Heaviside step function. For the Nernst 
signal we then get 



2CE,QAt)dt 



{Ey{t)QoiO)) dt. (62) 



Thus, in order to evaluate the response it is enough to 
consider the odd part of Q{t). 

Let us also make the following observation: At finite 
times the response will contain transients, which we will 
not be interested in and which do not contribute to the 
stationary response. Indeed, any contribution to Qoit) of 
the form of a total time derivative df /dt will contribute 
a term 

{Ey{0)^)dt = {Ey{0)f{^)) {EyiO)fi-^)) 

(63) 

to Cat, which vanishes provided f(t) is stationary in equi- 
librium. We will now treat the Langevin and RCSJ case 
separately. 



2. Langevin dynamics 

For the case of Langevin dynamics the dynamical ac- 
tion S obtains from substituting 77 using Eq. (|12p in the 
Gaussian noise probability distribution 

Ei:^^^'W^*) (64) 



resulting in 



^Lang. - / E {^^^ - ' (65) 
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where 



/i 86^ 2e 



and correspondingly 



(66) 



(67) 



Since = 2eVi/h is even and Ifj is odd under time re- 
versals^, the odd part of Q is 



where Ai — \j is the longitudinal part of the total current 
flowing through the links of the lattice, and ^ij = fii — fij 
is the transverse part, with the fij defined on the dual 
lattice sites, i.e., on the plaquettes, adjacent to the bond 
ij as in Fig. |4l The /i/:s are often referred to as loop 
currents. Neither of these contribute to the transport 
current As in the Langevin case we set /*J = 0, 
whereby the heat current equals the energy current. 

The dynamical action corresponding to Eqs. ((72H73)) 
can be expressed as 



54 



Qo{t) = ^^a;,M.F, 



El R 

ij) 



dt. 



1 h 



2T2 2e ^ ^ 

1 h 



2T^ 2e 



Y,M-x,9,)I^^. (68) 



Putting = {xi + Xj)/2 we may rewrite this as 



(74) 

The first term represents the Gaussian distribution of the 
white noise current J^" , and r]ij is a Lagrange multiplier 
to enforce the constraints Eq. ([7^ . In a temperature 
gradient the local temperatures are position dependent. 

The functional integration is over the variables 9, 7^" , 
rjij, Xi, and Integrating over I"j and rjij we get 



and 



The first term is a total time derivative of xf^ times a 
local energy e^j 

^ H ^^i = - ^ -^y cos 7,y- , (70) 

and will therefore not contribute to static response func- 
tions. The remaining part can be identified with the heat 
current in the Langevin model. 

This form of the heat current is directly formulated using 
the currents and potentials, and therefore simpler to use 
than Eq. (|52p . To show the equivalence of Eqs. ([71]) and 
([5^ it is necessary to go through some further steps. 
Before doing that, however, we consider the RCSJ case. 



3. RCSJ dynamics 

It is convenient to reformulate the equations of motions 
for RCSJ dynamics Eq. (IT9)) as 



(76) 

Again only the part which is odd under time reversal 
contributes to the heat current, 

2T'Qo = - ^ xt^V,, (X^ + It^ - A. + A, - , (77) 

{ij) 

since = Vij/Rij is even, while the other currents are 
odd. In this expression we may identify a contribution 



E ^"ij ^'^3 (-^^J + ^ij ) ~ E ^"iJ ' 



(ij) 



where 



2-K 

$0 



(78) 



(79) 



i^* = 4j-h/^+i^ + /U=i^ + A,-A,-HM.„ (72) 



is the local energy defined on the links of the lattice. 
Being a total time derivative this does not contribute in 
a stationary state. The remaining part can be rearranged 
into 

^^Ij^ij (-^^ - +Mii) 

^^^XiVi ^ (Ai - Aj + /ijj) 

i jeUi 



E 

jeMi 



Xi — A, — CoVi, 



(73) 



E 



^jj o(^ + ^0 ^ XijAtj 



{X^ - Xj + . 
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The second line equals Y^^XiCoViVi = -^Y^iXi^CoVi^, 
again a total time derivative which does not contribute. 
Finally, A.^ — Xj + = /*°' , so that the heat current for 
the RCSJ model becomes 



^s^7;(y^ + Vj)~x^^A,]ll°\ (80) 



4- Magnetization contribution 

Equations ^ and ([501) apparently differ. We 

will now show the equivalence of these formulations. We 
split the current into transverse and longitudinal parts 
as in Eq. (l46l) . with the /// defined on the dual lattice, 
whose sites are denoted by r/ = {xi,yi). With this it is 
possible to rewrite the last term in Eqs. (|7T|) . (|80|) . as 

~ $Z " ~ (^y ~ ^^i- 

(ij) 



The first term is a total time derivative of ^ 
where ef/ = i^y(Ai — A^) is the local magnetization 



M 



energy of the loops to ground [cf. Eq. (l49l) ]. The sec- 
ond term is a total time derivative of x^ef^ , involv- 
ing the magnetization energy ef^ — ^/i/^/ of the loops 
of the lattice. The remaining part corresponds exactly 
to the last term of Eq. ([52]). Note that for the mod- 
els discussed initially A = ($o/27r)A, i.e., only spatially 
uniform fluctuations are included in the vector potential. 
Then eff = ej = exactly, and Eqs. ([5^ and ([7T|). ([501) 
become identical (for open boundary conditions). In the 
more general case where local fluctuations are allowed 
they differ only by a total time derivative, which does 
not contribute to the transport coefhcients. 



C. Additional remarks 

The two derivations of the heat current given in 
Sec. IIII Al and IIIIBI above agree. The latter one shows, 
in addition, the validity of the standard Kubo formula 
Eq. for calculating the Nernst response of the trans- 
verse electric field to an applied temperature gradient (in 
absence of an electric transport current). The resulting 
expression should also give, via an Onsager relation, the 
response of the heat current to an applied transverse elec- 
tric current. This holds provided the electric transport 
current is subtracted from the current in Eqs. ([7T|). 
(|M)) . showing that Eq. (^51) is indeed the correct form. 

In both the Langevin and the RCSJ case the heat cur- 
rent is given by similar expressions, but in the RCSJ case 



the current /*°* includes also capacitative, resistive, and 
noise currents in addition to the supercurrent. 

As mentioned, Eq. ((25)) . being directly formulated in 
the currents and phases, have a clear advantage over 
(|52p. They are equivalent for systems with open bound- 
ary conditions along the temperature gradient. In sim- 
ulations periodic boundary conditions are convenient to 
use in order to eliminate surface effects. This, however, 
makes things more subtle as then the magnetization is 
not uniquely determined by the currents; Adding a con- 
stant to every fij does not change /*°* in Eq. (|^ . One 
logical possibility seems to be to impose an extra condi- 
tion on the average magnetization, e.g., define it to be 
zero at any moment, or fij — 0, so that no magne- 
tization contribution should be added in this case. An- 
other option is to fix one particular /i/, and in effect use 
Eq. ([25)) also for periodic boundary conditions. We opt 
for this latter condition, since it stays closer to the ex- 
perimental open system situation, while getting rid of 
surface effects. This choice can be further justified by 
comparing analytic results for open and periodic bound- 
ary conditions obtained in a spin wave approximation, to 
be discussed next. 



IV. HEAT CONDUCTIVITY AT LOW 
TEMPERATURE AND ZERO MAGNETIC FIELD 

At low enough temperatures and zero magnetic field 
the fluctuations will be small so that it is sufficient to 
consider linearized versions of the models introduced ear- 
lier. The only nonlinear circuit element is the Joseph- 
son junction, and by linearizing the Josephson relation 
I- J — Ifj sin jij w Ifjlij, the models are reduced to a 
network of capacitors, resistors and inductors, with ef- 
fective inductances Cij ~ {h/2elfj). In this spin wave 
approximation vortices are absent, hence there will be 
no Nernst effect, but the thermal conductivity will still 
be nonzero. For the Langevin case the resulting model 
can be mapped to one of heat conduction by phonons 
in a harmonic crystal coupled to local heat baths, which 
has an analytic solutioui^ For a 2D infinite square lat- 
tice the spin wave heat conductivity is independent of 
temperature and given by 



2ej 



sin^(7ra;) 



_ .1 _ 1 X fc^/c 
^ 4 it' 2e7 



4sin'(2f )+4sin^(^) 
kslc 



dxdy 



0.1817 



2e7 



(81) 



We now turn to the heat conductivity of the linearized 
RCSJ model on a 2D square array of size ^1 — Lx-^v 
and unit lattice constant. In this case, the heat current 
[Eq. §UI] includes the total current Iff = + + /f,- + 
/j", which is purely transverse due to Eqs. ([19| and ()20p . 
In this sum, I^j and //^ have no transverse component, 
so that /*°' = + /"-L. The transverse part of the 
supercurrent If^ is entirely due to vortices and vanishes 
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in the spin wave approximation, so that the total current 
only consists of the transverse component of the noise 
current, = More explicitly, this result can be 

derived from the equations of motion. On a square lattice 
it is convenient to label the links by the coordinate r and 
direction ii — x,y, and introduce forward and backward 
difference operators V^/(r) =V ^f{r -\- fi) — f{r + p,) ~ 
/(r). Introducing rescaled variables 0; = h9i/2e and 
A ~ hA/2e the equations of motion (|T9| . (|20l) are in this 
limit 



R 



(82) 



(83) 



Multiplying Eq. (|82|) with the lattice Green's function G 

(solving — VVGrr' = ^r,r') givCS 



(84) 



Using (IHll) and 
is 

tC 



the total current on a square lattice 



rtot ro I jH I fS I 



r/1 



- Vr, 



(85) 



We can here identify the longitudinal part of the noise 
current /r'J = - I]r'i/ ^M^rr' Vi./^'^, and the average 
(k = component) — TiJ2r ^rv^- The total current is 
thus just the transverse part of the noise current 

(86) 



rtot 



pi _ rn|| 
r/j, r/j, 



rnO 



The heat current in the x-direction can then be written 

as 

1 



(87) 



where = 2sin(fc^/2), = Y^^K^ = '^T^n^./L^, 
and = 0, . . . , L^ — 1. Since this is proportional to (5(i), 
the correlation function (x(i)x(O)) has to be evaluated at 
t = when inserted into the Kubo formula, i.e., it is given 
by an equilibrium correlation function. 

In the more general case the total current contains also 
the transverse part of the supercurrent, which is deter- 
mined by the vortices. Clearly the vortex contribution 
appears on top of the spin wave background calculated 
in this section. We may evaluate the correlation function 
{Xr^J.Xr'v) for the full model, including a capacitance Co 
to ground and cosine interaction. On a square lattice the 
RCSJ Hamihonian Eq. ^ is 



r r,/j, 

(91) 

where = 4>r and f{(j),A) is the "potential energy" 
involving the cosine interaction. Switching to Fourier 
space 

^ = ^ E l{Co+GK^)V^^V^+^CQiKl+kl)+fic^, A). 

(92) 

From here it is easy to calculate the required averages, 
since the partition function factorizes at the classical level 
and the averages are just Gaussian integrals. We get 



Co + CK^ ' 



Gn ''^ 



5^., (93) 



independent of /(</), A), so that 

(Xr.Xr'.> = ^ E e^'^-f"--'-') J |1 + e^"^ f (Vkl^^k) 

k 

+ (x + i)(x' + i)(AX), (94) 



Xr 



{x + -S,,)A,. (88) 



This result is intriguing because the dynamics of the sys- 
tem is completely independent of as seen from the 
equations of motion. Correspondingly, the correlation 
function which enters the Kubo formula Eq. ([23]) factor- 
izes: 

{jQit)jQio)) = l^ ^ {xr,m,^it)xMOK'tio)) 

= E Ur,{t)xr'.m{Ktm'i{Q))- (89) 

Furthermore, the transverse noise current correlation 
function 



2kBT , 
R 



X2 



(90) 



(XraXr'a) = XX' {AyAy I , 
{XrxXr'y) = {XryXr'x) = 0. 



(95) 
(96) 



Performing the sum over r, r' in ([5^ and using (|90p . 
the heat conductivity becomes 



,_ ks l^(l + K.^)(l-^ 



kT^O 



Co/C + 



RC 120. 



(97) 
(98) 
(99) 



The second term k" originates from (A^Aj^) and corre- 
sponds to the magnetization contribution. 

The resulting heat conductivity Eq. Wf}i has some no- 
table properties. Firstly, it is proportional to 1/RG so 
that it is well defined only for finite C. In this respect the 



12 



RSJ model, without shunting capacitors, i.e., with C — )■ 
is pathological. Secondly, when Cq = the sum over k 
in (|98|) is logarithmically divergent in the infinite system 
i — > oo (in 2D). For finite large L the heat conductivity 
has a logarithmic size dependence 



1 , L 

AttRC a ' 



(100) 



where a is the lattice spacing, k is thus not a bulk prop- 
erty of the RCSJ model in 2D. It is interesting to note 
that several other low dimensional models of heat con- 
duction display an anomalous size dependence, often tied 
to momentum conservationj2Si2I In the present case the 
diverging behavior is most likely due to the long range 
Coulomb interaction. A finite Co makes the charge- 
charge interaction exponentially small on distances larger 
than the screening length A = ay^C/Co, and also yields 
a system size independent Hgw when L ^ X. The calcu- 
lation above was done for periodic boundary conditions. 
We have repeated it for open boundary conditions in the 
x-direction. The difference is very small and tends to 
zero as l/L^ when increases. 

As discussed above, the vortex contribution appears 
on top of the temperature independent spin wave back- 
ground just calculated. Computationally it is often con- 
venient to project out the spin wave contribution by ex- 
cluding the transverse noise current /"-'- in Eq. (|5^ be- 
fore the averaging in Eq. (|23)) . 



V. SIMULATION METHODS 

The equations of motion for Langevin dynamics 
[Eq. ([T2|) and Eq. (fT3|) ] are solved numerically using a 
simple forward Euler discretization with a time step of 
At = 0.02. The RCSJ dynamical equations [Eq. (HH) and 
Eq. (|20p ] have a more complicated structure and are also 
second order in time, which makes the solution numer- 
ically more intensive. To solve these we use a leap-frog 
type discretization scheme, with time step At = 0.04. 
For a system of TV grains one then generally has to solve 
a system of A'^ -I- 1 coupled equations in each time step 
{N — 1 for the phases {9i} and 2 for the twists A). Note, 
however, that the equation system is sparse, so an effec- 
tive way to solve them is to employ an LU factorization 
algorithm, since the complexity of such a scheme goes as 
the number of nonzero entries, which are of the order of 
N here. This is far better than the direct method^^'^^ 
of multiplying with the lattice Green's function, which is 
here a dense matrix with ^ nonzero entries. In both 
the Langevin and the RCSJ case the sampling is per- 
formed during 4 • 10^ units of time, after an equilibration 
time of 10% of this. 

We have simulated systems of sizes up to 120 x 120, 
but except the case of k for RCSJ dynamics, finite size 
effects are unimportant for systems larger than L > 20, 
and thus only systems of size 20 x 20 are considered. 



The transport coefficients are calculated from the 
Kubo formulas Eqs. where the upper limit in 

the time integrals is replaced by a large enough time (> 
the correlation time), such that the cumulative value of 
the integral has stabilized its value. The Nernst signal e jv 
obtained from the Kubo formula in Eq. (|22p is also com- 
pared with the value of bat calculated in two other ways. 
The first is simply to apply a small temperature gradi- 
ent in the x-direction and measure the resulting electric 
field Ey , and the second is through the Onsager relation 
gn — Pxy/T, where /3 is the response of the heat current 
in the a:-direction to an applied electric current in the y- 
direction, as defined in Eq. ([5]). This gives an important 
consistency check of our calculations. A similar check is 
performed for the heat conductivity k. The Kubo for- 
mula turns out to be the most efficient way to calculate 
the response, since one does not have to worry about 
nonlinear effects. 



A. Units 

In the simulations we redefine time and temperature 
to make them dimensionless (we also measure length in 
units of a lattice constant a). For both Langevin and 
RCSJ dynamics, temperature is rescaled according to 



T ~^T- 



(101) 



The dimensionless time is obtained from the transforma- 
tion 



t t 



2ej' 



(102) 



for Langevin dynamics. In the RCSJ case the rescaling 
is 



t t- 



2eRr 



(103) 



(In the case where the junction parameters vary from link 
to link J"^, R and C denote a characteristic magnitude.) 
From this follows that for Langevin dynamics, the Nernst 
signal Cat is given in units ofks/^e^, the thermal conduc- 
tivity K in units of ksl'^ /2e^, and the resistivity p in units 
of ?i/(2e)^7. For RCSJ dynamics the Nernst signal tN is 
measured in units of 2ekBR/h, the thermal conductivity 
K in units of 2ekBRI'^ /h, and the resistivity p in units 
of the shunt resistance R. The dimensionless parame- 
ter = 2eR^I'^C /h (the ratio of the two times scales 
RC and h/2eRI^) controls the damping. For Q ^ 1 the 
system is underdamped and for Q ^ 1 it is overdamped. 



B. Time discretization of the heat current 

It is crucial to use a symmetric time discretization 
of the heat current, Eq. (1251) . in the numerics. For 
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Langevin dynamics, while it is sufScient to use a for- 
ward Euler discretization for the integration of the equa- 
tions of motion, the voltage Vi = h9i/2e appearing in 
Eq. (pSj) has to be approximated by a centered difference 
e{t) « {Oit + At) - e{t - At)}/2At. In the RCSJ case 
the total electric current is (after rescaling time) 

lT(f) = (sin7.,(t) + 7.,(i) + Q%,it)) + I^it). 

(104) 

In the symmetric leap-frog scheme we use, 9 is defined 
on integer time steps t = nAt, while the first order time 
derivative 9 is defined only on half-integer time steps t = 
(n -I- 1/2) At, so 9{t) has to be calculated as the average 
of 9 at the two adjacent time steps 

Ht) ~ ^{(){t + At/2) + 9{t - At/2)}. (105) 

The second order time derivative 9(t) is symmetrically 
defined as 

9{t) « ^{^(i + Ai/2) - 9{t - At/2)}. (106) 

The same applies for the twist variables A. These defini- 
tions make the heat current J^{t) as defined above nat- 
urally symmetric around t. An interesting aspect here is 
that by choosing the time step as At = 2Q'^, the RCSJ 
equations of motion discretized by the symmetric leap- 
frog scheme actually reduce exactly to the RSJ equations 
of motion discretized using an asymmetric forward Euler 
scheme. Moreover /*°*(i) becomes the sum of the super- 
, resistive, and noise currents discretized by a forward 
Euler scheme as it should for the RSJ model, while the 
voltages in Eq. (P5|) are kept symmetric due to the defini- 
tion Eq. (jlOSp . The RSJ model is therefore best thought 
of as a special case of an overdamped RCSJ model with 
= At/2 (with our choice of At = 0.04 this corre- 
sponds to = 0.02). 

We find that the heat current J^it) is very sensitive 
to the discretization used. In fact, it is critical to use 
the symmetric way of defining (t) to obtain consistent 
results when calculating the heat conductivity n either 
using a Kubo formula or by applying a small temperature 
gradient. 

VI. RESULTS AND DISCUSSION 

A. Zero field thermal conductivity 

Figure [S] shows the thermal conductivity k in zero 
magnetic field for fairly underdamped RCSJ dynamics 
{Q = 10) for different system sizes L. At low T it tends to 
the spin wave value given by Eq. (jlOOp . which in dimen- 
sionless units becomes k ~ {1 / inQ^) \n{L / a) . For large 
Q this background value is quite small. For smaller val- 
ues of Q the background increases and soon overwhelms 
the vortex contribution. In the following it will therefore 
be subtracted. 




T 



FIG. 5. (Color online) Heat conductivity k vs temperature T 
for &n LxL square lattice with underdamped RCSJ dynamics 
(Q = 10). The different curves correspond to different sys- 
tem sizes L. The inset shows the logarithmic dependence on 
system size L at low T (T = 0.1). The circles are simulation 
data and the smooth red curve is the analytic result obtained 
from a linearized model. 




T 



FIG. 6. (Color online) Heat conductivity k vs temperature 
T for an L X L square lattice with RSJ dynamics. Here 
the temperature independent spin wave background of n has 
been subtracted. The different curves correspond to differ- 
ent system sizes L. Inset: The same curves but scaled with 
1/ ln(L/0.7). The collapse is very good over the entire temper- 
ature range, except from close to the transition temperature 
Tbkt ^ 0.9. 



The dependence on system size L for low temperatures 
can be seen in the inset of Fig. [SJ The dependence is 
logarithmic and follows very well the form in Eq. (|100l) . 
shown as the red curve in the inset. In Fig.[B]we can see n 
as a function of temperature in the strongly overdamped 
limit C — >■ (corresponding to RSJ dynamics), but now 
with the harmonic spin wave background subtracted. In 
RSJ dynamics where C = 0, the spin wave contribution 
is formally proportional to (5(0), which translates to 1/ At 
in the numerics, where At is the time step used in the 
discretization (remember that RSJ dynamics with finite 
timestep At can be viewed as a special case of RCSJ 
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FIG. 7. (Color online) Nernst signal cn, heat conductivity 
K, and electrical resistivity p vs temperature T at different 
fillings / for a 20x20 square lattice with Langevin dynamics. 



dynamics with — At/ 2). What is left after subtract- 
ing this part can be interpreted as coming mainly ft'om 
the motion of vortices. The curves of n start out very 
small for low temperatures, but increase rapidly on ap- 
proaching the Berezinskii-Kosterlitz-Thoulesal^ temper- 
ature Tbkt — 0.9, where the unbinding of thermally in- 
duced vortex- antivortex pairs makes a large contribution 
to the thermal conductivity. At around T = 1.4 k reaches 
its maximum value, followed by a slow decrease for higher 
temperatures. Note that k, even after subtracting the 
background, shows a logarithmic dependence on the sys- 
tem size. The inset of Fig. [5] displays the curves for dif- 
ferent system sizes divided by the factor ln(i/a), with 
a = 0.7. The collapse of the curves onto a single one 
is very good over the entire temperature range (except 
from very close to Tbkt — 0.9, where small deviations 
are expectedly seen). 

In the spin wave approximation the logarithmic size de- 
pendence seems to be related to the long range Coulomb 
interation. Screening can be introduced by adding capac- 
itances to ground on every grain, which removes the di- 
vergent beha vior on length scales larger than the screen- 
ing length C/Cq. In a similar spirit, it seems likely 
that the logaritmic size dependence of the vortex contri- 
bution is tied to the unscreened Coulomb interaction in 
the RCSJ model. 




FIG. 8. (Color online) Nernst signal ejv, heat conductivity 
K, and electrical resistivity p vs temperature T at different 
fillings / for a 20x20 square lattice with RCSJ dynamics {Q 
— 0.5). In the k plot, the temperature independent spin wave 
background has been subtracted. 

For Langevin dynamics, on the other hand, k is only 
weakly dependent on system size and quickly converge to 
a size independent bulk value (for L > 20). Moreover, 
the finite size effects of e^r and p are also negligible for 
L > 20, for both Langevin and RCSJ dynamics. We will 
therefore stick to lattices of size 20 x 20 in the remainder 
of this paper. 

At low T for Langevin dynamics (see the lowest curve 
in the middle panel of Fig. [7]) k goes to the spin wave 
value Eq. (|81l) , but decreases slightly upon increasing the 
temperature until temperature induced vortices become 
plentiful near the BKT transition, where k increases 
again reaching a maximum around T sa 1.25 and then 
starts to decrease. 



B. Low fields 

We now turn to the case of a relatively weak applied 
transverse magnetic field. In Fig. [7] and H] we see a col- 
lection of simulation results for fillings / = to 0.05 on 
a 20 X 20 square lattice for Langevin and RCSJ dynam- 
ics (Q = 0.5) respectively. The filhng / = B (Api) /$o, 
where (Api) is the average plaquette area, represents 
the average number of magnetic field induced vortices 
per plaquette in the system. The lowest nonzero filling 



15 



/ = 0.0025 = 1/20^ corresponds to one field induced vor- 
tex in the system. Increasing the filhng, i.e., raising the 
magnetic field will cause the vortex density to increase, 
but as long as the typical vortex separation is much larger 
than the coherence length ^ can be thought of as the 
short distance cut-off or lattice spacing a in our model) 
the effects of discreteness are negligible and the model 
should describe a continuous two-dimensional (or quasi- 
two-dimensional) type-II superconductor. The case of 
strong magnetic fields are discussed in Section IVI CI 

Focusing on the Nernst signal cat in Fig. [7] and [HI we 
notice a very steep increase of cat at low temperatures 
to a maximum between T = 0.10 and 0.15. At higher 
temperatures slowly decreases and the tail persists 
up to about T — 2, which is roughly twice the BKT 
transition temperature Tbkt ~ 0.9. These features are 
qualitatively similar for both Langevin and RCSJ dy- 
namics. The main difference between the models in the 
shape of Cat is the plateau-like part of the curves present 
in the RCSJ case around T = 0.7 for low fillings. The 
peak height increases rapidly with filling, before it starts 
to decrease for higher fillings. The position of the Nernst 
signal peak depends slightly on the filling / and moves 
towards higher temperatures for larger fillings. 

The qualitative features of the Nernst signal are in 
agreement with experiments on several superconduc- 
torsJ-i^i^i^ A detailed comparison is, however, possi- 
ble only if the temperature and magnetic field depen- 
dences of the parameters 7, i?, /"^ are taken into ac- 
count. For Josephson-j unction arrays made up of tun- 
nel junctions the Ambegaokar-Baratoff formula'*'^ ~ 
(7rA(T)/2ei?)tanh(A(T)/2fcBT) can be used. In bulk 
superconductors is proportional to absolute square of 
the superconducting order parameter jV'P, which goes as 
~ (TMF _ j.-)^ ^^^^ j.MF^ r^Yie situation in the high-T^ 
cuprates is more complicated, since a model for the re- 
laxation rate, 7(T)^^ in the Langevin case, or i?(T) in 
the RCSJ case, is also required. The thermoelectric co- 
efficient axy = &N I P may have an advantage here, since 
both Gat and p are proportional to the relaxation rate, 
which therefore drops out.— iS. We will return to a^y in 
the next subsection. 

In the second row of Fig. [7] and [8] the heat conductivity 
K is plotted as a function of temperature at different fill- 
ings / (in the RCSJ case the spin wave background has 
been subtracted) . Note first how similar the low temper- 
ature part of the curves of At are to the Nernst signal. 
The onset and the peak positions of the two quantities 
agree to a high degree. For Langevin dynamics (Fig. [7]) 
K is finite in the limit T — > 0. For the two lowest fillings 
/ = 0.0025 and 0.005, k increases quickly as a function 
of T and then falls off slowly below the T — value of 
~ 0.20 to suddenly increase again around Tbkt and reach 
a second maximum followed by slow a decrease at higher 
temperatures. At fillings above / = 0.005 the thermal 
conductivity follows the same pattern, but does not fall 
below the T — > value until temperatures above Tbkt- 
In the high temperature regime T > 1.5 all curves, re- 



gardless of filling /, fall onto a single curve. The curves 
for overdamped RCSJ dynamics with Q — 0.5 (Fig. 
share this feature of two maxima. In this case, how- 
ever, the falloff after the first maximum is even more 
pronounced and persists up to higher fillings, at least 
/ = 0.05. Note that the temperature independent back- 
ground contribution to k has been subtracted in this fig- 
ure. 

The double-peak behavior seems to indicate two sep- 
arate contributions to the thermal conductivity at high 
and low temperatures. The first maximum at low T is 
probably caused by the increased mobility of the field in- 
duced vortices, which also gives the sharp rise of cat. This 
contribution diminishes with increasing T, until the un- 
binding of temperature induced vortex- antivortex pairs 
around Tbkt — 0.9 makes k large again. This latter con- 
tribution totally dominates the previous one at higher 
temperatures, causing curves for different fillings to con- 
verge. 

Looking at the resistivity p, it displays, not the same 
but a qualitatively similar T dependence, for Langevin 
and RCSJ dynamics. For Langevin dynamics the effect 
of varying filling / is somewhat more apparent, and the 
rise at Tbkt is also a bit steeper than for RCSJ dynamics. 

The results for the RCSJ model discussed previously 
have been for the overdamped case Q — 0.5. Upon re- 
ducing the damping and moving into the underdamped 
regime, cat is effectively unchanged in the low tem- 
perature region (apart from a trivial change of scale). 
However, k and p change slightly at high temperatures 
T > 1.5, in that k decays somewhat faster and the falloff 
starts at a lower temperature, while p increases more 
quickly as function of T, than in the underdamped case. 



1. Vortex heat transport 

In our models there is no Hall effect {axy = 0), lead- 
ing to bn = axy I <J XX = OixyPxx, and so we can obtain yet 
another transport coefficient, the off-diagonal component 
of the thermoelectric tensor Uxy from e^v and p. As men- 
tioned one advantage with axy is that, unlike ejv and k, it 
does not depend on the time constants 7 and R^^^ This 
makes comparison with experimental data easier. Fur- 
thermore, from phenomenological theories of vortex mo- 
tion axy has the interesting interpretation of the entropy 
per vortex)^ suggesting that axyT can be identified with 
the transported heat per vortex. This can also be seen 
from the Onsager relation dxy = axyT together with the 
definition of the transverse electrothermal conductivity 
axy = Jx /E^y > which is measured under conditions with 
no applied temperature gradient. Now, since Ey is pro- 
portional to the transverse vortex current through the 
relation Ey — ^oJx, the quantity J§ /Ey is just the ratio 
between the heat current and the vortex current, i.e., a 
measure of the average transported heat per vortex. 

Another way of calculating the heat transported per 
vortex, or more precisely the ratio between the heat cur- 
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FIG. 9. (Color online) A plot of k/cn and a^yT vs temper- 
ature T for RCSJ dynamics {Q = 0.5) on a square 20 x 20 
lattice at filling / = 0.02. Both of these measure the heat 
transport per vortex (J^/JJf$o). Inset: The / dependence 
of these two quantities at low T (T = 0.1). The two smooth 
green curves are fits to the form a + & In / (light green: a = 
0.11, h = -0.19; dark green: a = -0.14, b = -0.22). 




ted as functions of temperature T. While they are not 
equal, they do agree very well at low T (and for low fill- 
ings /), so here one can approximately speak about trans- 
ported heat per vortex. n/eM is consistently larger than 
a,j;yT, and the difference grows with increasing T (and in- 
creasing /). A clue to why this happens can be found con- 
sidering the difference in what drives the vortex motion 
in the two cases, as mentioned above. When applying a 
temperature gradient, heat can be transported even with- 
out a net flow of vortices, being mediated solely through 
the interactions between vortices at different tempera- 
tures. This is not possible when the driving force is the 
transverse electric current, and naturally explains why 
k/cjv is always greater than a^yT. The magnetic field 
(or filling /) dependence of k/bn and a^yT at a low 
fixed temperature is shown in the inset of Fig. [HI The 
plot is in lin-log scale and clearly shows a logarithmic de- 
pendence at low fillings, which is quite similar for both 
quantities. Such a logarithmic dependence obtains from 
an ideal gas treatment of the vortices, the Sackur- Tetrode 
entropy per vortex being ~ — In /. The vortices are, how- 
ever, strongly interacting. A crude way to estimate the 
interaction effects on the transport entropy would be to 
assume that the available volume per vortex is reduced 
by a factor of N, the number of vortices. This then gives 
a contribution Inn/N — — In / to the configurational en- 
tropy, i.e., also a logarithmic dependence. 

Also notice that in the low temperature region of 
Fig. [HI up to about T ~ 0.5, k/cn and a^yT are only 
weakly temperature dependent. This means that axy 
roughly falls of as ^ 1/T for low temperatures. A fit to 
- 1/T of a^y for RCSJ dynamics (Q = 0.5) at / = 0.01 
is displayed in Fig. [TUl The inset in log-log scale reveals 
that for temperatures above T ~ 0.8, axy falls off much 
faster, somewhere close to ~ 1/T*. These features are 
valid also for Langevin and RSJ dynamics, as well as for 
other types of lattices. 



FIG. 10. (Color online) The off-diagonal component of the 
thermoelectric tensor axy — eiv/p vs temperature T for RCSJ 
dynamics (Q — 0.5) on a square 20 x 20 lattice at filling 
/ = 0.01. At low T the curve follows the power law ~ 1/T. 
In the high temperature region the falloff is much faster, about 



rent and the vortex current, is to simply divide the vor- 
tex contribution k of the longitudinal heat conductivity 
(i.e., with the spin wave background subtracted) with the 

Nernst signal, k/cn — ( _v t )/( -v''t ) ~ l^v Note, 
however, that in this context the ratio / Ey measures 
the transported heat per vortex in a system driven by 
a temperature gradient Va;T, as opposed to the case of 
o-xyT = J§ /Ey, where the driving force is a transverse 
electric current Jy. An equality of k/cat and a^yT would 
imply a strict proportionality of the heat current on 
the vortex current J^. 

In Fig.|n]K/ejv and a^yT from our simulations are plot- 



C. Intermediate and high fields - effects of 
granularity 

Going to higher magnetic fields the type of lattice 
structure starts to play an important role, as geomet- 
ric frustration will affect vortex transport in this regime. 
The models at hand are then more valid as descriptions 
of granular superconductors. 



1. Square lattices 

Figures [TT] and [T^ show simulation results as a function 
of filling / for Langevin and RSJ dynamics (corresponds 
to RCSJ dynamics with Q = y/om) on a 20 X 20 square 
lattice. At low fillings the Nernst signal e n and the heat 
conductivity k both show a sharp increase culminating 
in a maximum around / = 0.05 — 0.15, depending on 
temperature, followed by a decrease up to half-filling. 
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FIG. 11. (Color online) Nernst signal ejv, heat conductivity k, 
and electrical resistivity p vs filling / at different temperatures 
for a 20x20 square lattice with Langevin dynamics. 
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FIG. 12. (Color online) Nernst signal ejv, heat conductivity k, 
and electrical resistivity p vs filling / at different temperatures 
for a 20x20 square lattice with RSJ dynamics. 



This "tilted-hill" profile of the Nernst signal seems to 
be generic and is found in a number of experiments on 
cuprates and ordinary type-II superconductorsj^^ Note 
how cn is always zero at / = and 1/2, due to vortex - 
vacancy symmetry. The heat conductivity k on the other 
hand stays finite at / = 1/2. On a perfectly periodic lat- 
tice all physical quantities are, because of the symmetry 
of the XY model Hamiltonian [Eq. (fTTj) ]. periodic in fill- 
ing /, with period one, and also mirror symmetric (for p, 
k) or antisymmetric (e^v) around / = 1/2. Thus, all in- 
formation is contained in the region / = — >■ 0.5, which 
is displayed here. 

Now, lowering the temperature the curves have signif- 
icant structure due to geometric frustration as the filling 
is varied through different commensurate values. At fill- 
ings such as / = 1/8, 1/5, 1/3, 2/5, all three transport 
coefficients cat, k and p are reduced due to vortex pinning 
to the underlying lattice. This is particularly apparent 
at the 2nd lowest temperature T ~ 0.2 (cyan colored 
curves) at f = 1/3 and 2/5. The observant reader may 
also have noticed that the Nernst signal actually goes 
negative in a region below half-filling. In fact a small 
region of negative Nernst signal appear also right below 
/ = 1/3, and it is plausible that this occurs below other 
commensurate fillings as well, over certain temperature 
intervals. This sign reversal of cat is a new effect^^ seen in 
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FIG. 13. (Color online) Contour plot of the Nernst signal ejv 
for a 20x20 lattice with Langevin dynamics in the tempera- 
ture - filling plane. The dashed white line joins the maxima 
of ejv. 



all our simulations independent of the type of dynamics 
used (Langevin or ovcr-/underdamped RCSJ) and per- 
sists also for lattices with moderate geometric disorder 
(see e.g. Fig. [15]). As discussed in a previously published 
paper of ours,^^ a negative vortex Nernst signal [given 
the definition in Eq. ^] implies vortex transport in the 
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FIG. 14. (Color online) Nernst signal ejv, heat conductivity k 
and electrical resistivity p vs filling / at T = 1 for random lat- 
tices of size 20 X 20 with dynm = 0.8, using Langevin dynamics. 
Results shown are for two models with different critical cur- 
rent distributions Jf^ — I'^ (red) and /fj ~ (black). Each 
curve is an average over eight disorder realizations. 



direction opposite of heat transport. We argue that this 
is possible, since around these special fillings there is a 
temperature regime, where mobile vortex vacancies can 
exist on top of a pinned vortex lattice. The vortex vacan- 
cies then diffuse down the applied temperature gradient, 
creating a net vortex flow in the opposite direction and 
thus a negative Nernst signal. 

Comparing the Nernst signal versus / for Langevin and 
RSJ dynamics in Fig. [Til and fT2l respectively, one sees an 
almost exact agreement (as opposed to bn versus T at 
low fillings, where Langevin and RSJ dynamics are less 
similar). Also increasing the damping parameter Q for 
RCSJ dynamics, i.e., going from the overdamped to the 
underdamped Hmit results in hardly any changes in the 
qualitative filling dependence of cat , k and p. This should 
indicate that these quantities in this regime are governed 
by geometric frustration effects and are rather insensitive 
to model specific details. 

As a summary of the Nernst effect on a square lattice 
(for Langevin dynamics) we provide in Fig. [T3] a contour 
plot of ejv in the temperature - filling plane. The red 
regions indicate a large Nernst signal and the blue ones 
a signal which is close to zero, or even negative (the dark 
blue blob in the upper left corner). 
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FIG. 15. (Color online) Nernst signal ejv, heat conductivity k, 
and electrical resistivity p vs filling / at T = 1 for random lat- 
tices of size 20 X 20 with different degrees of disorder set by the 
parameter dmin using Langevin dynamics. Each curve is an 
average over 16 disorder realizations. The lower panel shows 
examples of Voronoi lattices obtained from random packings 
of grains with different dmin. 



2. Random lattices 

In order to model random granular superconductors we 
have carried out simulations on randomly connected net- 
works as defined in Sec. Ill CI In Fig. [T3] we compare e^r, 
K and p as functions of filling at T = 1 for two different 
models with Langevin dynamics on a random lattice with 
dmin = 0.8. In the first model (red curve) we set the criti- 
cal currents of every junction to a constant Ifj = I'^. The 
second model (black curve) has the critical currents pro- 
portional to the contact area (or strictly contact length 
djj in 2D) between the grains, Ifj ~ djj, see Fig. [31 

In a geometrically disordered system without perfect 
periodicity the Nernst signal and the resistivity are no 
longer periodic as a function of filling. We have there- 
fore extended the curves up to / = B {Ap{) /^q = 2 
{{Ap\) =1/2 for the random lattices). The Nernst signal 
shows an even steeper increase at low / and peaks even 
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FIG. 16. (Color online) Contour plot of the Nernst signal ejv 
in the filling - grain size standard deviation plane for random 
lattices with Langevin dynamics at T = 1. 

earlier than in the square lattice case. We also see that 
Cat is somewhat reduced in the model with If, ~ dj-4 . The 
negative region is narrower but above / = 1 the curves 
are essentially the same. Looking at the heat conduc- 
tivity the two models display quite similar behavior as 
a function of temperature, although in the second model 
the amplitude of k is effectively halved. A more dramatic 
difference can be seen in p though, where the model with 
added disorder of the critical currents seems almost like 
a smoothed-out version of the first one, for which p has 
more dramatic features at high fillings. At very high 
fillings / ^ 1 the flux through each plaquette modulo 
<I>o becomes approximately random in [0, $0] making the 
Nernst signal vanish and k and p constant. The results 
for RCSJ dynamics show a qualitatively similar behavior, 
with some minor quantitative differences. 

Results for different strengths of randomness (different 
values of dj^in, see Sec. Ill C| are shown in Fig. [151 From 
this figure it is also apparent that much of the structure 
in e^v, K and p as a function of filling is reduced as we 
move from less towards more disordered lattices. (The 
apparent lack of visible geometric frustration effects at 
specific fillings are, however, in the case of these random 
lattices due to the fact that we average over many (8 or 
16) different disorder realizations.) 

Fig. [TC] summarizes the Nernst effect in our random 
lattices. It displays a contour map of the Nernst signal 
in the filling - grain size standard deviation plane. The 
large blue region at the bottom corresponds to a negative 
Cat, whereas the red region to the left represents a large 
positive Nernst signal. 

VII. CONCLUSIONS 

We have modeled heat and charge transport in two- 
dimensional granular superconductors. We have consid- 



ered regular square arrays and randomly connected net- 
works. Square arrays can be artificially fabricated using 
lithography, while random arrays of varying degrees of 
disorder may occur naturally in inhomogeneous thin-film 
superconductors. We consider two models for the dy- 
namics of the superconductors, relaxational Langevin dy- 
namics and RCSJ dynamics. An expression for the heat 
current is derived from these models. For the Langevin 
dynamics the heat current expression (|25)) is consistent 
with previous expressions derived microscopically or from 
time-dependent Ginzburg-Landau theory j^i^i^^ For the 
RCSJ model, however, it differs in that it contains the to- 
tal current, including the capacitative and resistive cur- 
rents, which shunt the supercurrent through the junc- 
tion. The expressions are used in numerical simulations 
to calculate the heat conductivity, the resistivity, and the 
Nernst signal. We find an anomalous logarithmic size de- 
pendence in the heat conductivity for the RCSJ model in 
zero magnetic field. This type of dependence is present 
also in the spin wave approximation valid at low T. We 
also find k to be divergent in the limit when the shunt- 
ing capacitor goes to zero, showing that the RSJ model 
without capacitors is pathological from the point of view 
of heat conduction. The Nernst signal and the resistivity 
are still well behaved in this limit. 

From our numerical simulations, we further find a 
highly nontrivial nonmonotonous temperature depen- 
dence in K at low magnetic fields in both models. In this 
regime granularity appears to have a negligible influence 
on the transport properties, and our results should apply 
also to two-dimensional phase-fluctuating bulk supercon- 
ductors. For low temperature T -C Tbkt ~ 0.9 and mag- 
netic fields it is possible to define the transported heat 
per vortex, which is found to depend logarithmically on 
filling /, while being approximately temperature inde- 
pendent. Note, however, that in our phase-only models 
the vortices are coreless. At higher T thermally excited 
vortices and antivortices start to infiunce the results and 
dominate the response. 

At higher fields, granularity becomes important and 
geometric frustration strongly influences the Nernst sig- 
nal, heat conductivity, and resistivity, leading to a highly 
intricate magnetic field dependence as shown in Figs. 1111 
These signatures should be possible to obtain di- 
rectly in experiments on regular Josephson- junction ar- 
raySj_similar to those carried out for the resistivity in 
Ref. |43| . or in patterned thin film superconductors. 
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